function [bestever, array, total_time] = saea_rfs(func_num, m, ts_num, K, sub_size, nIter, maxfe)
% SAEA-RFS 
% 基于代理模型的进化算法 (SAEA) + 随机特征选择 (RFS)
% 分解高维问题为多个子问题，并利用差分进化 (DE) 
% 子问题。通过代理模型减少真实函数评估次数。
%
% INPUTS:
%   func_num - 目标函数编号 (目前支持 1 或 15)
%   m        - 种群大小
%   ts_num   - 初始训练样本数量
%   K        - 子问题数量
%   sub_size - 子问题最大维度
%   nIter    - 子问题优化最大迭代次数
%   maxfe    - 最大函数评估次数
%
% OUTPUTS
%   bestever   - 全局最优解 (目标函数最小值)
%   array      - 优化历史记录 [累计评估次数, 当前最优值]
%   total_time - 总运行时间 (目前默认 0)


% 屏蔽 RBF 矩阵可能出现的奇异值警告
warning('off', 'MATLAB:nearlySingularMatrix');
warning('off', 'MATLAB:singularMatrix');

global initial_flag
initial_flag = 0; % 用于 cec2013 测试函数初始化标志

%% 问题维度设置
if func_num > 12 && func_num < 15
    D = 905; % 特殊重叠函数的维度
else
    D = 1000; % 默认维度
end

%% 搜索空间边界设置
if ismember(func_num, [1,4,7,8,11,12,13,14,15])
    lu = [-100*ones(1,D); 100*ones(1,D)];
elseif ismember(func_num, [2,5,9])
    lu = [-5*ones(1,D); 5*ones(1,D)];
elseif ismember(func_num, [3,6,10])
    lu = [-32*ones(1,D); 32*ones(1,D)];
end
d = D;  % 问题维度

%% 初始化档案 archive (训练样本)
XRRmin = repmat(lu(1, :), ts_num, 1);
XRRmax = repmat(lu(2, :), ts_num, 1);
archive = XRRmin + (XRRmax - XRRmin) .* lhsdesign(ts_num, d); % 拉丁超立方采样

% 计算真实函数值，CEC2013期望列向量格式
archive_fitness = cec2013(archive(:,1:d)', func_num);
archive(:, d+1) = archive_fitness';

FES = ts_num; % 当前函数评估次数
cycle = 0;
flag = 'cubic'; % RBF 类型，后续可改为 TPS、Gaussian 等

bestever = min(archive(:,d+1)); % 初始全局最优
fprintf('cycle = %d\tBest fitness: %e\n', cycle, bestever);
array = [FES, bestever]; % 历史记录初始化

%% 初始化种群 Parent
XRRmin = repmat(lu(1, :), m, 1);
XRRmax = repmat(lu(2, :), m, 1);
Parent = XRRmin + (XRRmax - XRRmin) .* lhsdesign(m, d); % 拉丁超立方生成初始种群

%% 主优化循环
while (FES < maxfe)
    
    %% 1. 使用 RFS 构建子问题并训练代理模型
    ts = archive; % 当前训练样本
    [lambda, gamma, dv, sp] = rfs(sub_size, ts(:, 1:d), ts(:, end), flag, K);
    % lambda/gamma: 子问题的 RBF 模型参数
    % dv: 子问题决策变量索引
    % sp: 子问题训练样本索引
    
    %% 2. 子问题优化
    [bestY, index] = min(archive(:, end)); % 找到全局最优的训练样本
    bestX = archive(index, 1:d); % 全局最优解
    sub_eval_min = []; % 用于存储每个子问题的最优预测值及解
    
    for i = 1:K
        % 提取子问题种群
        sub_Parent = Parent(:, dv{i});   
        % 利用 RBF 模型预测子问题种群适应度
        fit_Parent = rbf_eval(sub_Parent, ts(sp{i}, dv{i}), lambda{i}, gamma{i}, flag);
        % 子问题搜索空间边界
        subBound = lu(:, dv{i});    
        
        % 使用 DE 优化该子问题
        for iter = 1:nIter
            % DE 生成新候选子代
            sub_child = de(sub_Parent, bestX(1, dv{i}), subBound);
            
            % 利用 RBF 模型预测子代适应度
            fit_child = rbf_eval(sub_child, ts(sp{i}, dv{i}), lambda{i}, gamma{i}, flag);
            
            % 合并父代与子代并排序
            fit_pop = [fit_Parent; fit_child];
            [~, y] = sort(fit_pop); % 升序排序
            com_pop = [sub_Parent; sub_child];
            
            % 环境选择：选择前 m 个最优个体
            sub_Parent = com_pop(y(1:m), :);
            fit_Parent = fit_pop(y(1:m));
        end
        
        % 保存子问题最优解及其预测值
        sub_eval_min(i, 1) = fit_Parent(1);                 % 最优适应度
        sub_eval_min(i, 2:length(dv{i})+1) = sub_Parent(1, :); % 最优解
        
        % 更新整体种群对应子问题维度
        Parent(:, dv{i}) = sub_Parent;  
    end
    
    %% 3. 确定使用真实函数评估的解
    solution = bestX; % 先使用全局最优作为初始解
    [min_model, sign] = min(sub_eval_min(:, 1)); % 找到预测最优的子问题
    temp = rbf_eval(Parent(1, dv{sign}), ts(sp{sign}, dv{sign}), lambda{sign}, gamma{sign}, flag);
    
    if temp <= min_model
        % 若代理模型预测优于历史最优，则采用 Parent 中的子问题解
        solution(:, dv{sign}) = Parent(1, dv{sign});
    else
        % 否则使用子问题的预测最优解
        solution(:, dv{sign}) = sub_eval_min(sign, 2:length(dv{sign})+1);
    end
    
    %% 4. 搜索新的解并更新档案
    diff = setdiff(solution, archive(:, 1:d), 'rows'); % 新解排除重复
    if ~isempty(diff)
        % 真实函数评估新解
        diff_fitness = cec2013(diff(:, 1:d)', func_num);
        diff(:, d+1) = diff_fitness';
        FES = FES + size(diff, 1); % 更新函数评估计数
        archive = [archive; diff]; % 更新档案
    end
    
    %% 5. 更新历史记录和全局最优
    cycle = cycle + 1;
    bestever = min(min(archive(:, d+1)), bestever);  % 更新全局最优
    array = [array; FES, bestever];                 % 记录历史
    fprintf('cycle = %d\tBest fitness: %e\n', cycle, bestever);
end

total_time = 0; % 暂不统计时间，可根据需要添加 tic/toc

% 恢复警告设置
warning('on', 'MATLAB:nearlySingularMatrix');
warning('on', 'MATLAB:singularMatrix');

end % function
